%% Assignment 8.2
% CIE4440

close all;clearvars;clc;

%% 8.2.0 Initiate data

load data82;    % loading Q and time

startDate   = [1994,9,1,0,0,0];

data82.timeseries   = timeseries(data82.Q);
data82.timeseries.Name = 'Discharge Rhine at Lobith';
data82.timeseries.TimeInfo.Units = data82.timeUnit;
data82.timeseries.TimeInfo.StartDate = datestr(startDate);
data82.timeseries.DataInfo.Units = data82.QUnit;

figure(1)
plot(data82.timeseries,':b');

%% 8.2.1 Use the Chezy formula to estimate the stage and velocity at every hour
% Define constants
Ch      = 30;
ib      = 1/1000;
width   = 250;

h   = (data82.Q./(Ch*width*sqrt(ib))).^(1/1.5);

timeseries   = timeseries(h);
timeseries.Name = 'estimated stage Rhine at Lobith';
timeseries.TimeInfo.Units = data82.timeUnit;
timeseries.TimeInfo.StartDate = datestr(startDate);
timeseries.DataInfo.Units = 'h';

figure(2)
plot(timeseries,':b');

figure(3)
cax     = gca;
set(cax,'NextPlot','add');
title('Not taking hysteresis into account');
xlabel('Q m3/s');
ylabel('h m');
plot(data82.Q,h);

%% 8.2.2 Assume c=3/2u
u       = data82.Q./(width*h);
c       = 1.5 * u;

%% 8.2.3 Calculate dh/dt and 1/ci
dh      = [0;diff(h)];  % Forward differencing
dt      = 3600;
dhdt    = dh./dt;

ci_1    = (c*ib).^-1;

%% 8.2.4 Motivate your answer on whether or not hysteresis is important to consider in this case for this river.
% You can for instance assume that the stage that you calculate is the true stage. Then calculate the discharge Q with the Jones equation and plot this Q against h. 

% I reccon it is important to include hysteresis in your calculations,
% since the bedshape changes for different stages. The floodplanes form a
% large factor.

Qjones  = data82.Q.*sqrt(1+ci_1.*dhdt);

figure(4)
cax     = gca;
set(cax,'NextPlot','add');
title('Taking hysteresis into account');
xlabel('Q m3/s');
ylabel('h m');
plot(Qjones,h);

warning('no difference!')

